OPEN 3 ACCESS Freely available online 



•0-PLOS I ONE 



Increasing fMRI Sampling Rate Improves Granger 
Causality Estimates 

Fa-Hsuan Lin 1 ' 2 *, Jyrki Ahveninen 3 , Tommi Raij 3 , Thomas Witzel 3 , Ying-Hua Chu 1 , liro P. Jaaskelainen 2 , 
Kevin Wen-Kai Tsai 1 ' 2 , Wen-Jui Kuo 4 , John W. Belliveau 3 

1 Institute of Biomedical Engineering, National Taiwan University, Taipei, Taiwan, 2 Department of Biomedical Engineering and Computational Science, Aalto University, 
Espoo, Finland, 3 Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Boston, Massachusetts, United States of 
America, 4 Institute of Neuroscience, National Yang-Ming University, Taipei, Taiwan 



CrossMark 



Abstract 

Estimation of causal interactions between brain areas is necessary for elucidating large-scale functional brain networks 
underlying behavior and cognition. Granger causality analysis of time series data can quantitatively estimate directional 
information flow between brain regions. Here, we show that such estimates are significantly improved when the temporal 
sampling rate of functional magnetic resonance imaging (fMRI) is increased 20-fold. Specifically, healthy volunteers 
performed a simple visuomotor task during blood oxygenation level dependent (BOLD) contrast based whole-head inverse 
imaging (Inl). Granger causality analysis based on raw Inl BOLD data sampled at 100-ms resolution detected the expected 
causal relations, whereas when the data were downsampled to the temporal resolution of 2 s typically used in echo-planar 
fMRI, the causality could not be detected. An additional control analysis, in which we SINC interpolated additional data 
points to the downsampled time series at 0.1 -s intervals, confirmed that the improvements achieved with the real Inl data 
were not explainable by the increased time-series length alone. We therefore conclude that the high-temporal resolution of 
Inl improves the Granger causality connectivity analysis of the human brain. 
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Introduction 

Determining causal mechanisms by which different brain areas 
interact to support cognition and behavior has been a persistent 
challenge in neuroscience. Whereas analyzing synchronization of 
cerebral activations can identify cortical areas acting in concert, 
revealing causal influences among them requires measures of 
effective connectivity [1-3]. Previously, effective connectivity analyses 
of human PET [4] and fMRI [5-8] data have been conducted 
using structural equation modeling (SEM), which aims at 
determining directional modulations across the activated areas 
using covariance or correlation matrices derived from the 
measured time series [9]. However, a major limitation of SEM 
is that it requires strong a priori assumptions on the number and 
directionality of connections, which are often difficult to justify or 
validate. Similar limitations exist in dynamic causal modeling 
(DCM), which also requires a priori models of directional 
connections [10-12]. To circumvent such limitations, the tech- 
nique of Granger causality [13] has been applied to data obtained 
with both EEG [14-22] and fMRI [18,23-29]. The main 
advantage of Granger causality over SEM and DCM is that it 
can estimate the directionality of modulations within a network 
without a priori assumptions on which connections are active and 
on directions of the connections. Essentially, Granger causality 



tests how additional information improves prediction of the future 
of a given time series. In other words, a Granger causal influence 
from a time series "X" to time-series "Y" exists if the combined 
information from "X" and "Y" predicts the future of "Y" better 
than information from "Y" alone. 

Functional MRI of the human brain [30] with blood 
oxygenation level dependent (BOLD) contrast [31,32] is the 
prevailing method for studying brain functions noninvasively. 
There are two major limitations to using BOLD fMRI for 
causality modelling. First, BOLD signals are vascular responses 
that lag the underlying neuronal events by seconds [33] and may 
show notable voxel-to-voxel latency variability at the individual 
level [34]. However, it has been suggested that with appropriate 
modelling to obtain neuronal activity estimates, BOLD fMRI can 
still be used for causality modelling [35] . The other challenge for 
using BOLD fMRI in Granger causality estimation is the rather 
low sampling rate, which is critically important in all time series 
modeling. Typically fMRI Granger causality analyses use data 
sampled at the rate of approximately 1-2 s [24,26-29]. Such a 
slow sampling rate, which is necessary to achieve whole-brain 
fMRI coverage at a spatial resolution of 3-4 mm, provides only 
about 10-15 samples during the 20-30 sec duration of a canonical 
hemodynamic response function [36]. Estimating Granger cau- 
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sality from fMRI time series recorded at such a low sampling rate 
can be problematic. 

Using the recently developed dynamic functional magnetic 
resonance inverse imaging (Inl), one can achieve an order of 
magnitude faster sampling rate. Inl is based on the utilization of 
simultaneous measurements from multiple channels of an RF head 
coil array, and by solving sets of inverse problems Inl can detect 
dynamic changes of the BOLD fMRI signals at ~10 Hz sampling 
rate with whole-brain coverage and approximately 5-mm spatial 
resolution at the cortex [37-39]. Our recent study suggests that, 
Inl hemodynamic responses can elucidate neuronally related 
timing information when cross-subject and within-region variabil- 
ity is suppressed by averaging [40]. 

Several studies have consistently suggested that the sensitivity 
and stability of Granger causality values can be critically improved 
if the temporal sampling rate is high enough [26,41-47]. However, 
to our knowledge, there have been to date no studies empirically 
demonstrating this. Based on our data showing that the BOLD 
fMRI signal can reflect neuronal timing at the group level [40], 
here we hypothesize that increasing the fMRI sampling rate using 
Inl one can provide more robust and sensitive Granger causality 
estimates compared to conventional multi-slice EPI acquisitions. 
We test this empirically using Inl fMRI with 10-Hz Inl sampling 
rate and a simple visuomotor detection task, which generates feed- 
forward inter-area information flow [48]. Three different time 
series were used in this study: fMRI raw time series, hemodynamic 
response function after General Linear Model, and the estimated 
neuronal activity using hemodynamic deconvolution. Time series 
with lower sampling rates (2 Hz, 1 Hz, 0.5 Hz, and 0.2 Hz) were 
artificially generated by either discarding Inl samples or interpo- 
lating sub-sampled time series in order to keep the same number of 
time points. Our results indicate that the high sampling rate 
provided by Inl can robustly improve detection of causal 
modulations between cortical areas. 

Materials and Methods 

Ethic statement, subjects, and tasks 

This study was approved by the Institute Review Boards of 
National Taiwan University Hospital and National Yang Ming 
University. Written informed consent approved by the Institute 
Review Board of National Taiwan University Hospital and 
National Yang Ming University was obtained from each subject 
prior to participation. Healthy human volunteers (re = 23, 6 
females, all right-handed, age 22-30 years) were presented with 
left or right visual hemifield reversing (8 Hz) checkerboard stimuli 
in a rapid event-related fMRI design. The hemifield checkerboard 
subtended a 4.3° visual angle and was generated from 24 evenly 
distributed radial wedges and eight concentric rings of equal 
width. The stimuli were presented using the Psychtoolbox [49,50]. 
Stimulus duration was 500 ms; the onset of each presentation was 
randomized with a uniform distribution of inter-stimulus intervals 
varying from 3 to 16 s (average 10s). Part of data has been 
previously reported in studying the correlation between hemody- 
namic and neuronal activity [51]. 

The subjects were instructed to press the button upon detecting 
a visual stimulus, presented randomly at the left or right side of the 
screen, with the hand ipsilateral to the stimulus. Thus, there were 
two conditions: right visual hemifield-right hand (R condition) and 
left visual hemifield-left hand (L condition). This relatively simple 
task was chosen for its feedforward connectivity from visual to 
motor systems. Accordingly, our a priori hypothesis predicted 
directional information flow from visual to motor cortices. 



Twenty-four stimulation epochs were presented during four 
240 s runs, resulting in a total of 96 trials per subject. 

Structural MRI acquisitions and reconstructions 

Structural T x -weighted MRIs were acquired with a 3T scanner 
(Tim Trio, Siemens, Erlangen, Germany) and a 32-channel head 
phased array coil using a standard MPRAGE sequence (repetition 
time/echo time/inversion time [TR/TE/TI] = 2,530/3.49/ 
1100 ms, flip angle = 7°, partition thickness = 1.33 mm, image 
matrix = 256x256, 128 partitions, field-of-view = 2 1 cmx21 cm). 
The location of the gray-white matter boundary for each 
participant was estimated with an automatic segmentation 
algorithm to yield a triangulated mesh model with approximately 
340,000 vertices [52-54]. This cortical model was then used to 
facilitate mapping of the structural image from native anatomical 
space to a standard cortical surface space [52,53]. Between- 
subjects averaging was done by morphing individual data through 
a spherical coordinate system [55]. 

fMRI inverse imaging (Inl) acquisitions and 
reconstructions 

BOLD-contrast fMRI data were acquired using inverse imaging 
[39,56], which included a reference scan to collect spatial 
information from different channels in a radio-frequency coil 
array and a set of dynamic scans to achieve high temporal 
sampling rate (>10 Hz) with whole brain coverage. The Inl 
reference scan was collected using a single-slice echo-planar 
imaging (EPI) readout, after exciting one thick coronal slab 
covering the entire brain (FOV 256 mmx256 mmx256 mm; 
64 x64 x64 image matrix) with the flip angle set to the gray matter 
Ernst angle of 30° (considering the Tl of gray matter is 1 second 
at 3T). Partition phase encoding was used to obtain spatial 
information along the anterior-posterior axis (Y direction). EPI 
readout had frequency encoding along the superior-inferior (Z 
direction) and phase encoding along the left-right axis (X 
direction). We used TR= 100 ms, TE = 30ms, band- 
width =2604 Hz and a 12.8-s total acquisition time for the 
reference scan, allowing the coverage of the whole-brain volume 
with 64 partitions and two repetitions. 

For the Inl functional scans, we used the same volume 
prescription, TR, TE, flip angle, and bandwidth as for the Inl 
reference scan. The principal difference was that, to achieve the 
high temporal resolution, partition phase encoding (in the Z 
direction) was removed so that the full volume was excited, and the 
spins were spatially encoded by a single-slice EPI trajectory, 
resulting in a coronal X/Z projection image with spatially 
collapsed projection along the anterior-posterior direction. The 
i-space Inl reconstruction algorithm [57] was then used to 
estimate the spatial information along the anterior-posterior axis. 
In each run, we collected 2,400 measurements after collecting 32 
measurements in order to reach the longitudinal magnetization 
equilibrium. A total of 4 runs of data were acquired from each 
participant. The total fMRI acquisition time for each subject was 
16 minutes. 

Inl data were reconstructed time-point by time-point using the 
minimum-norm estimate method [39,56], which generated 2,400 
reconstructed volumes in each run. Subsequently, we used 
General Linear Model (GLM) to identify functional areas. The 
design matrix of GLM included the impulse trains of stimulus 
onset convolving the finite-impulse-response (FIR) basis function, 
which included 6-s pre-stimulus baseline and 24-s post-stimulus 
response, as well as DC and linear drift terms to estimate the 
hemodynamic responses. For statistical analyses, the noise levels in 
the reconstructed images were estimated by calculating the 
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standard deviation of the time series during the 6-s pre-stimulus 
interval. Using these noise estimates, dynamic statistical paramet- 
ric maps (dSPMs) were derived as the time-point by time-point 
ratio between the Inl reconstruction values and the baseline noise 
estimates. Such dSPMs can be assumed to be t distributed under 
the null hypothesis of no hemodynamic response [58]. 

Each subject's Inl time series, including the estimated hemo- 
dynamic responses and raw fMRI time series before GLM, were 
spatially registered into their cortical surface space by using a 12- 
parameter affme transformation between the volumetric Inl 
reference scan and the MPRAGE anatomical space (FSL, 
http://www.fmrib.ox.ac.uk/fsl). The resulting spatial transforma- 
tion was subsequently applied to each time point of the 
reconstructed Inl volume. For inter-subject averaging the individ- 
ual functional results were co-registered to a spherical brain 
surface representation [59]. 

Regions-of-interest (ROIs) were determined by the spatial 
distribution of the temporally average t statistics greater than 4.0 
(Bonferroni corrected />-value<0.05) between 4 s and 7 s after 
visual stimulus onset. This interval was chosen to capture the 
maximum BOLD response. Time courses for each ROI were 
extracted and averaged within the ROI for each subject 
separately. 

Echo-planar imaging (EPI) collection and analysis 

We also collected EPI with 2 s TR in order to use empirical 
data to compare the detection sensitivity to causality modulation 
using Inl data. The subjects were instructed to perform the same 
lateralized visuomotor task as that during the Inl acquisitions. 
Data were collected from 13 subjects over one 4-minute run, 
where 30 left hemifield and 30 right hemifield visual stimulations 
were randomly presented. EPI parameters were: TR/TE = 2000/ 
30 ms, field-of-view (FOV) = 220x220 mm, matrix = 64x64, slice 
thickness = 4 mm, flip angle = 90°. For each subject, thirty-four 
trans-axial slices with no gap were acquired with the spatial 
coverage of cerebrum and cerebellum. 

EPI fMRI data were first pre-processed for motion correction, 
slice timing correction, and spatial smoothing (12 mm 3D full- 
width-half-maximum Gaussian filter) by using the FreeSufer 
software (http://surfer.nmr.mgh.harvard.edu). We used GLM to 
identify visual and motor cortices with the FIR bases described 
above. Estimated hemodynamic activity from each subject was 
spatially registered to a common surface coordinate system using 
FreeSufer. The BOLD signals were then averaged across subjects 
and the dynamic significance of the BOLD signal was calculated 
using the dSPM procedure described above. 

Time series preparation 

In this study, we analyzed the causal modulation using the fMRI 
time series directiy, instead of using the estimated HRF. In 
preparation of these fMRI time series, GLM was first used to 
identify the location of visual and sensorimotor cortices (see section 
fMRI inverse imaging (Inl) acquisitions and reconstruc- 
tions and Echo-planar imaging (EPI) collection and 
analysis above). These time series were then averaged within 
each ROI from each subject. The DC value and the linear drift 
were also removed from these time series. 

This study evaluates how the increased temporal resolution in 
Inl acquisitions and reconstructions can help elucidate the causal 
interactions in the human visuomotor system. This was done by 
parametrically downsampling the Inl data to generate time series 
from 0.1 -s to 0.5, 1, and 2-s temporal resolution, the last of which 
is quite common in whole-head fMRI studies. In addition, to avoid 
the confound of different sampling rates and thus different lengths 



of time series, we also generated SINC -interpolated time series 
from the sub-sampled time series such that all time series have the 
same length in the causality analysis. Finally, we also used EPI 
time series (before GLM) to estimate causality modulations in 
order to compare the results based on down-sampled Inl data. 

Granger causality and conditional Granger causality 
analysis 

We used an auto-regressive (AR) model to implement the 
Granger causality modeling. Consider a zero-mean time series at 
the "destination node", y(t). The^ th -order AR model description of 

y(fi is: 

p 

y(t)= ^2a k y(t-k)+s y (t), t=l---n, (1) 

k=\ 

where are AR model coefficients, and Ey(t) is the residual time 
series of AR model fitting at the destination node, n is the total 
number of samples in the time series. Provided with the time series 
from a "source", x(t), we can model the bivariate time series of x[t) 
and _)>(/) as: 

\y{t\m T = 

p (2) 
Y,A k \y(t-k),x(t-k)] T + [e y , x (t),e XrX (t)] ,t=l - n 

k = l 

where A k 's are the AR model coefficient matrices, and \e y JJ), 
£ x,x(^)] T is the joint bivariate residual time series of AR model 
fitting at the destination node and source node. The Granger 
causality metric is 

GC x ^ y = log s 2 y (t) I g ejjt) j (3) 

Since the bivariate time series \y(t), *(/)] T contains the information 
of univariate time series y(t), from Cauchy inequality we can 
conclude that G My is well defined and positive since the quotient 
inside the logarithm is greater or equal to one. Previously, it has 
been suggested that the inference of G My can be calculated by 
referring ((n—p)/p) (exp(G My )-(n— 2p)/(n— p)) to an F distribution 
with (p, n—2p) degrees of freedom [60]. 

For each pair of ROI time series ("X" and "Y"), we respectively 
calculated the Granger causality G My and G y ^ x . Previous 
simulation studies suggest that naive computation of Granger 
causality over fMRI signals can be misleading [26] . However, the 
influence difference (the difference between the pair G x _^ y and 
G y ^ x ) can be much more robust [25,26]. Following this rationale, 
for each pair of ROIs, we only calculated the difference between 
two Granger causality values in order to identify the dominant 
influence direction. 

The implementation of the AR modeling was done by the 
ARFIT algorithm [61,62], where the optimal model order was 
jointly determined by model fitting (i.e., favoring the model with 
smaller power of the residual time series after fitting) and model 
parsimoniousness (i.e., favoring a lower order model). This 
implementation has been previously used to explore the causal 
modulation of epileptic spike propagation measured by MEG [63]. 
Allowing the AR model ordering ranging between 1 and 20 or the 
largest order allowed by the number of sample, we used the 
Bayesian information criterion to choose the optimal order of the 
AR model. In the Granger causality analysis on time series of a 
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Table 1. Timing indices and the full-width-half-maximum (FHWM) of the group-average hemodynamic responses in five regions- 
of-interest. 









V 


PCC 


PreM 


S 


M 


Onset (s) 


L 


0.63 


0.91 


0.69 


1.04 


1.18 




R 


0.60 


1.03 


0.84 


1.04 


1.19 


Time-to-half (s) 


L 


1.70 


1.90 


1.90 


2.20 


2.40 




R 


1.70 


2.00 


2.20 


2.20 


2.40 


Time-to-peak (s) 


L 


3.60 


3.70 


4.00 


4.10 


4.50 




R 


3.70 


4.10 


4.20 


4.20 


4.40 


FWHM (s) 


L 


4.60 


4.30 


5.10 


4.60 


4.70 




R 


4.80 


4.80 


4.70 


4.60 


4.60 



V: visual cortex; PCC: parietal cortex; PreM: pre-motor cortex; S: sensorimotor cortex; M: motor cortex. L: left hemisphere (R condition). R: right hemisphere (L condition). 
doi:1 0.1 371 /joumal.pone.01 0031 9.t001 



pair of ROIs, we used the minimum of the estimated AR model 
orders to estimate the GC for the sake of model simplicity. 

We used a non-parametric approach to estimate the statistical 
significance of the influence difference, because the null distribu- 
tion of such influence difference has no analytic form. Further- 
more, if we want to use the F/ Chi-square distributions in causality 
estimation, samples in the time series need to be independent. This 
may not be the case in BOLD-contrast fMRI time series with 
different sampling rates. Our non-parametric analysis started from 
generating the null distribution of influence difference using the 
Adjusted Amplitude Fourier Transform (AAFT), a method of 
generating surrogate time series while preserving the linear 
correlation structure of the original time series and the marginal 
distribution [64]. The permutation procedure was repeated 1000 
times for Granger causality analysis and 100 times for conditional 
Granger causality analysis, because the latter analysis took a longer 
computation time. We defined the p-value as the number of 
occurrences of the Granger causality values using a swapped 
source time series exceeding the Granger causality value using the 
original source time series. The p-value of the Granger causality 



difference in the group analysis was calculated accordingly as the 
ratio between the number of the Granger causality difference 
values from permuted time series higher than that from the 
original time series in each subject and the total number of 
permutation test across subjects. Note that such a ^-value 
calculation amounts to a fixed-effect analysis. 

Granger causality can be complicated by a "common source" 
problem. Specifically, the causal modulation between a pair of 
time series can be mediated through other time series within the 
set of analyzed time series. This problem has been partially 
addressed by the conditional Granger causality approach [25,65- 
67] . In the case of studying the causal modulations between the 
time series pair \y(t), x(t)] T , suppose that a multivariate auto- 
regressive model can be used to describe all other time series 
[^(Z)] T . The conditional Granger causality cG x _^ y is calculated as 
the ratio between the residual variance in the joint time series \y(t), 
z(t)] T and the residual variance in the joint time series [x(C), y{tj, 




Figure 1. (Left) Locations of the five ROIs (f statistics of the BOLD signal averaged between 4.0 s and 7.0 s after the visual stimulus 
onset) in each hemisphere. (Right) Hemodynamic time courses and estimated neuronal activity using hemodynamic deconvolution at five ROIs. 
doi:1 0.1 371 /journal.pone.01 0031 9.g001 
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Table 2. The optimal order of the AR model of the time series at five ROI's. 
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PCC 


PreM 
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TR = 0.1 s 
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11 


11 


11 


11 


12 




R 


10 


8 


9 


9 


9 


TR = 0.5 s 


L 


3 


3 


3 


3 


3 




R 


4 


3 


4 


4 


4 


TR = 1 s 


L 


1 


1 


2 


1 


1 




R 


2 


2 


2 


2 


2 


TR = 2 s 


L 


1 


1 


1 


1 


1 


R 1 1 1 11 



V: visual cortex; PCC: parietal cortex; PreM: pre-motor cortex; S: sensorimotor cortex; M: motor cortex. L: left hemisphere (R condition}. R: right hemisphere (L condition). 
doi:1 0.1 371 /joumal.pone.01 0031 9.t002 



\y(t)At)\ T = 

£ B k \y(t-k),z(t-kf+ [e y , z (0,£zM \t=\---n 
[x{t\y(t),z(tf = 

£ C k [x(t - k),y(t ~ k)Mf ~kf+ [£*,,z W^wzM^wW] r , 

k=\ 

t=\---n 

cGC x ^ y = log ( £ 4(f) / £ ) 



Heuristically, conditional Granger causality evaluates the improve- 
ment of the additionally explained variance in the target time series 
by adding the source time series after controlling the information in 
all other time series. Statistical significance can be similarly derived 
from the empirical null distribution. 

All calculations were done using Matlab (Mathworks, Natick, 
MA, USA). 

Results 

The visuomotor task elicited strongest BOLD signal at the left 
hemisphere, contralateral to the visual stimulus /motor response. 
Thus, results for the L condition are reported in the right 
hemisphere and the R condition in the left hemisphere. We 
identified five regions-of-interest (ROIs), including visual cortex 
(V), parietal cortex (PPC), pre-motor cortex (PreM), somatosen- 
sory cortex (S), and motor cortex (M) for both L and R conditions, 
as shown in Figure 1 . Notably, the details of the average time 
courses (Figure 1) show sequential BOLD activity with different 
onsets, time-to-half (TTH), time-to-peak (TTP) timing, and width 
of regional hemodynamic response (Table 1). 

Figure 2. shows that increasing the fMRI sampling rate 
significantly improved the sensitivity of our Granger causality 
estimates, as indicated by the emergence of multiple directional 
influences consistent with the visuomotor task. The dominant 
direction of information flow between any two ROIs (denoted X 
and Y) was determined by calculating the difference (Gx-^y - 



Gy-h>x) between the two uni-directional Granger causality 
estimates [23,26]. Table 2 lists the order of the optimal AR 
model for each ROI time series. P-values of all directions of 
information flow between ROFs were listed in Table 3. At the 
highest temporal resolution (TR = 0. 1 s), significant bottom-up 
causal influences were observed from visual to PPC, premotor, 
somatosensory, and motor cortices in both left and right 
hemispheres. Reducing the temporal resolution clearly decreased 
the number of significant causality influences. At TR = 0.5 s, left 
hemisphere still shows three strong feedforward connections. 
However, the significance levels of two connections (from visual 
cortex to premotor and motor cortices) have decreased. Further 
lowering the sampling rate down to 1 s and 2 s, the conventional 
TR for whole-head EPI, shows no significant feedforward 
connections. In addition to the influence differences, we also 
report the Granger causality values and their associated signifi- 
cance in Table 4, showing that almost all causality estimates were 
significant at TR< = 1 s and non-significant at TR = 2 s and that 
the majority of statistically significant GC connections were 
bidirectional. 

Granger estimates can be confounded by one source driving 
multiple targets. To reduce confounds related to potential 
common sources, we also calculated the difference conditional 
Granger causality values among pairs of the time series in 5 ROIs. 
Different from the Granger causality, conditional Granger 
causality estimates the specific information flow between two 
chosen ROIs while the information from all other ROIs through 
direct or indirect information flow are removed [65,67,68]. P- 
values are marginally different between Granger causality analysis 
and conditional Granger causality analysis in all sampling rates 
(Table 5). Again, increasing the TR gradually reduced the 
number of significant paths - no significant causal modulation was 
observed in both hemispheres at TR = 1 s and 2 s. 

To confirm that the increased sensitivity in detecting feedfor- 
ward connections is not due to different time series length, we 
SINC interpolated the time series from the subsampled data such 
that all time series had the same length before Granger causality 
analysis. For example, TR = 0.5 s data were first down-sampled by 
5-fold from the densely sampled time series with TR = 0. 1 s and 
subsequently SINC interpolated by 5-fold. Results in Figure 3 
shows that, while interpolation can artificially help improve the 
detection, slow sampling rates at the order of 1 s cannot provide 
significant feedforward connection estimates in left and right 
hemispheres. Specifically, we found that the significance levels for 
the V-^M and V^-PreM connections in the left hemisphere 
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Table 3. P-values of all directions of information flow estimated by Granger causality between ROI's at different sampling rates 
from 1000 bootstrap iterations. 
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Each cell contains two p-values. The left p-value is for the left hemisphere (R condition) and the right p-value is for the right hemisphere (L condition). Only lower- 
triangular off-diagonal entries are listed. Upper-triangular off-diagonal entries can be calculated by (1 -listed p-value). V: visual cortex; PCC: parietal cortex; PreM: pre- 
motor cortex; S: sensorimotor cortex; M: motor cortex. L: left hemisphere (R condition). R: right hemisphere (L condition). 
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decreased (increased /(-value) at TR = 0.5s. A significant 
V-^PreM connection shows up at TR = 0.5 s in the right 
hemisphere after data interpolation. Results with TR = 1 s 
suggested possible V-^PreM and V— »M connections. However, 
the trend of losing detecting feedforward connections at a slower 
sampling rate prevails. In particularly, TR = 2 s, the typical whole- 
head EPI sampling rate, shows no significant connections in both 
hemispheres. 

Finally, in addition to using down-sampled Inl data, we also 
used EPI data with TR = 2 s record the subjects's BOLD signal at 
visual and motor cortices in the same lateralized visuomotor task. 
Consistent with the down-sampled Inl simulations (Table 3), we 
found that there was no significant V— >M causal modulation in 
both R and L conditions (/(-values = 0.27 and 0.22 for R and L 
conditions, respectively). 



Discussion 

Our results demonstrate that the inverse imaging method with a 
20-fold increase in fMRI temporal resolution can gready improve 
the detection power of causality analysis in the human brain. 
Using whole-head Inl acquisitions with 100 ms temporal resolu- 
tion, we were able to infer directional causal influences between 
five functional areas activated during a lateralized visuomotor task 
in each hemisphere with both Granger causality and conditional 
Granger causality analyses. Consistent with the a priori predicted 
pattern of functional connectivity in a visuomotor choice-reaction 
task [48], the results were dominated by feedforward influences 
from visual to sensorimotor, premotor, and posterior parietal 
regions. These observations strongly underline the importance of 
high temporal sampling rates in determining effective connectivity, 
because the connectivity patterns only emerged at the 1 00 ms TR 



PLOS ONE | www.plosone.org 



6 



June 2014 | Volume 9 | Issue 6 | e100319 



Fast fMRI Improves Causality Estimates 



V V V V 



V V V V 



V V V 



u 
a. 
a. 



u 
a. 
a. 



o 



V V V V 



V V V V 



% 
o 

a. 



o 

II 



O r- O t- 



o\ i- 



% 



u 
a. 
a. 



vo t— 
o r\i 
»- o o 



U 
a. 
a. 



V V V V 



V V V V 



r\l ro m rM 



o 



o 

a. 



IN 

II 



Oi <N CO vfj 
fN r- ^- 



VD 


LA 


LD 


o> 




ro 






rM 


o 


d 


d 






cn 




LA 


ex* 




ro 




O 


d 


o 



rM t- t- 



*- rs <- 



PLOS ONE | www.plosone.org 



7 



June 2014 | Volume 9 | Issue 6 | e100319 



.en hemisphere Kiyru nemispiieie 

(R condition) (L condition) 



Fast fMRI Improves Causality Estimates 



R = 0.1 s 




p < 0.05 
0.05 < p < 0.1 



R = 0.5 S 




TR = 1 s 




TR = 2 S 




Figure 2. The dominant information flow calculated from the difference between two uni-directional Granger estimates among the 
visual (V), PPC, premotor (PreM), somatosensory (S), and motor (M) cortex ROIs at TR = 0.1 s, 0.5 s, 1 s, and 2 s. Significant causal 
modulations (p<0.05) were shown in thick yellow arrows and connections showing a trend of causal modulation (0.05<p<0.1) were shown in thin 
yellow arrows. 

doi:1 0.1 371 /journal.pone.01 0031 9.g002 



(Figures 2 and 3). Controlling the number of time points by 
interpolating the time series with a lower sampling rate can only 
partially mitigate the problem of losing sensitivity in correct 
causality estimates (Figure 3). 

The advantage of higher temporal resolution has been suggested 
by simulation studies [26,41-46], but has not been previously 
shown with empirical data. Moreover, the previous simulation 
studies assumed that regular EPI would be used, where any 
increase in fMRI temporal resolution has to be traded off for 
smaller spatial coverage, poorer spatial resolution, or jittered 
designs that greatly increase the duration of the imaging session 
[2,69]. Specifically, limited by the gradient slew rate, the 
acquisition of each echo-planar imaging slice takes approximately 
80 ms at ~3 mmx3 mm spatial resolution. Assuming that 30 
slices would be needed to cover the entire brain, a jittered stimulus 



design would require a 30-fold increase in data acquisition time, 
leading to impractical session durations (several hours). With the 
Inl approach in the present study, we obtained a whole-head 
volume in 100 ms, but it is possible to speed up the acquisition 
even further. The temporal resolution of Inl is determined by TE 
(30 ms) optimized for the BOLD contrast and desired field-of- 
view/ spatial resolution. Using partial Fourier acquisition and a 
higher readout bandwidth, the readout time could be reduced 
from 32 ms (2 KHz bandwidth and 64 lines) to 16 ms (6/8 partial 
Fourier, 48 lines at 3 KHz bandwidth) at the cost of reduced SNR. 
However, if using an echo-shifting pulse sequence [56,70], it is 
possible to achieve 20 ms TR at 30 ms TE. These future 
developments could allow for testing causal modulations between 
hemodynamic time series with putative latency differences well 
below 100 ms. 
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Table 5. P-values of all directions of information flow estimated by conditional Granger 
sampling rates from 100 bootstrap iterations. 
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Each cell contains two p-values. The left p-value is for the left hemisphere (R condition) and the right p-value is for the right hemisphere (L condition). Only lower- 
triangular off-diagonal entries are listed. Upper-triangular off-diagonal entries can be calculated by (1 -listed p-value). V: visual cortex; PCC: parietal cortex; PreM: pre- 
motor cortex; S: sensorimotor cortex; M: motor cortex. L: left hemisphere (R condition). R: right hemisphere (L condition). 
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In this study we chose Inl for its whole-brain coverage, 
reasonable spatial resolution at cortex, and high sampling rate 
(TR = 0.1 s). Note that recently there are other fast fMRI 
acquisition methods, such as generalized Inl [71], fast fMRI 
[72], MR-encephalography (MREG) [73], and simultaneous 
multi-slice EPI [74] methods. Each method has different 
magnetization excitation strategies and image reconstruction 
algorithms. Yet all of them can provide a sub-second sampling 
rate. As we demonstrate that fast fMRI acquisitions can help 
improve detecting causal modulations in human brain in general, 
researchers can choose any method adaptively under different 
theoretical or practical concerns without losing the detection 
power. 

Inl trades off a small amount of spatial resolution in one 
encoding direction for a substantially higher temporal resolution. 



Depending on the reconstruction methods, Inl has approximately 
5 mm spatial resolution at cortex using a 32-channel head coil 
array at 3T [37-39], which is in the range of spatial filtering/ 
smoothing that is typically applied in echo-planar imaging as a 
post-processing step. Without reaching the limit based on the 
electromagnetic theory [75,76], using a head coil array of more 
channels at a higher field (for example, 7T) could further improve 
the spatial resolution of Inl. Novel reconstruction methods 
targeted at suppressing the point spread function, such as using 
the minimum L-l norm constraint [77], can also be used to 
improve the spatial resolution in studies where the highest possible 
spatial resolution is critical. 

In addition to the neurophysiologically expected feedforward 
modulation from V-^PPC^preM^-M, our Inl results also 
suggest direct feed-forward connectivity from V to preM, M, 



PLOS ONE | www.plosone.org 



9 



June 2014 | Volume 9 | Issue 6 | e100319 



Fast fMRI Improves Causality Estimates 




Figure 3. The dominant information flow calculated from the difference between two uni-directional Granger estimates among the 
visual (V), PPC, premotor (PreM), somatosensory (S), and motor (M) cortex ROIs at TR = 0.1 s, 0.5 s, 1 s, and 2 s. The time series at 0.5 s, 
1 s, and 2 s were SINC interpolated such that all time series are of the same length as the 0.1 s time series, which contains the real measured data. 
Significant causal modulations (p<0.05) were shown in thick yellow arrows and connections showing a trend of causal modulation (0.05<p<0.1) 
were shown in thin yellow arrows. Importantly, only very little improvement is observed in the estimates where the number of observations is 
artificially increased using the SINC interpolation procedure: the strongest connectivity patterns are, clearly, observed in the 0.1 -s condition with only 
real data. 

doi:1 0.1 371 /journal.pone.01 0031 9.g003 



and S. Although these influences can be supported by animal 
models showing direct structural connections from visual to 
sensory and motor cortices in rats [78], it is also possible that 
different neurovascular coupling in the different ROIs caused 
smearing of neuronal temporal information that reduced the 
sensitivity of GC to detect the dominant feedforward connections 
in our estimates. 

In contrast to the present finding, previous studies have shown 
that EPI data may show causal modulation in the human 
sensorimotor system using EPI with TR in the order of 2 s 
[24,26]. However, this seeming discrepancy may be explained by 
the different nature of the present two-choice reaction-time task 



and those used in many previous studies. Here, the group average 
reaction time was less than 400 ms [40]. Therefore, without 
considering the differential vascular responses at visual and 
sensorimotor cortices, the BOLD signal time series in the visual 
and sensorimotor cortices were expected to be delayed by only a 
few hundreds of milliseconds. Such a latency is difficult to be 
resolved by regular EPI with TR = 2 s. Consequently, one can 
expect that quite similar BOLD signal time series in visual and 
sensorimotor cortices were elicited during this task. Two similar 
time series are actually difficult to use Granger causality analysis to 
reveal any causal modulation between them. In a simple 
theoretical example with two identical time series, it is clear that 
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Figure 4. The EPI (left column) and Inl (right column) time series at the visual and sensorimotor cortices from a representative 
subject (top panel) and the residual EPI (left column) and Inl (right column) time series at the sensorimotor cortex after AR 
modeling using the sensorimotor cortex time series alone and both sensorimotor and visual cortices (bottom panel). 
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the residual time series in modeling one time series will not be 
further improved by providing the other time series, because all 
information has been provided by the first time series already. In 
such a case, the Granger causality value will be insignificant 
because no further reduction on the power of residual time series 
after modeling the first one. 

To support our argument above, Figure 4 shows an example of 
two highly correlating (r— 0.446) EPI time series at visual and 
sensorimotor cortices from one representative subject. Figure 4 
also shows the residual time series after modeling the sensorimotor 
time series when either only the sensorimotor cortex time series 
was provided, or both the sensorimotor and visual cortices time 
series were provided. The residual time series power only changed 
marginally (108.0^-105.6; 2.3% reduction), in line with our main 
result of no significant Granger causality modulation between 
visual and sensorimotor cortices as measured by EPI with 
TR = 2 s. For comparison, we also show the Inl time series in 
visual and sensorimotor cortices from one representative subject, 
and the residual time series after modeling the sensorimotor time 
series when either only the sensorimotor cortex time series was 
provided, or both the sensorimotor and visual cortices time series 
were provided. Visually, it is difficult to discern the precedence of 
either time series. However, numerically the variance of residual 
time series at the sensorimotor cortex was apparently reduced 
when the visual cortex time series was provided (0.0026^0.0023; 
1 1.5% reduction). 

Finally, it should be noted that BOLD-contrast fMRI measures 
the vascular responses secondary to neuronal events. On top of 
information reflecting neuronal activity, there are well document- 
ed vascular confounds in fMRI time series [34,79]. Our previous 



study has shown that, with group averaging and improved 
sampling rate, inter-regional hemodynamic responses can be 
significantly correlated with inter-regional neuronal activity [51]. 
In this study, we chose a relatively simple two-choice reaction-time 
task with a priori assumption of observing strong feed-forward 
connectivity from the visual to the motor cortices. Such findings 
are consistent with a previous study [80] suggesting that fast 
sampling can eliminate confounding effects of differential HRF 
delays over areas by fine features of the HRF waveform (see also 
Table 1). However, while our results suggest that it is feasible to 
increase the fMRI sampling rate to enhance sensitive of causality 
estimates, caution must be always exercised when interpreting the 
results provided by these methods in the context of tasks eliciting 
more complex feed-forward/ feedback information flow patterns. 

Taken together, our results suggest that using MR Inl with a 
100 ms sampling interval (20-fold faster than conventional EPI), 
the sensitivity of detecting causal connectivity is significandy 
improved. We expect that this method can be used in other fMRI 
experiments to reveal effective connectivity when the vascular 
confound of the BOLD-contrast fMRI is carefully controlled and 
thus potentially open up entirely new possibilities for non-invasive 
imaging of effective connectivity in the human brain. 
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